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Abstract In this paper we analyze a 1+3 dimensional solution of relativistic hydrodynamics. We calculate 
momentum distribution and other observables from the solution and compare them to measurements from 
the Relativistic Heavy Ion Collider (RHIC). We find that the solution we analyze is compatible with the 
data. In the last several years many numerical models were tested, but it is the first time that an exact, 
parametric, 1+3 dimensional relativistic solution is compared to data. 
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1 Introduction 

\ In the last several years it has been revealed that the mat- 

■ ter produced in the collisions of the Relativistic Heavy Ion 
Collider (RHIC) is a nearly perfect fluid Q], i.e. it can be 
described with perfect fluid hydrodynamics. There was a 
long search for exact hydrodynamic solutions (solutions of 

, the partial differential equations of hydrodynamics) and 

■ several solutions proved to be applicable. There are 1+3 
dimensional solutions, as well as relativistic solutions - 
but no 1+3 dimensional and relativistic exact solution 

' has been tested yet. In this paper we extract observables 
from the relativistic, ellipsoidally symmetric solution of 
ref. [5] . We calculate momentum distribution, elliptic flow 
and correlation radii and compare them to RHIC data. 



2 Perfect fluid hydrodynamics 

Perfect fluid hydrodynamics is based on local conservation 
of entropy or number density (n) , energy- momentum den- 
sity (T^), expressed by so-called conservation equations: 

d^nu") = 0, (1) 
d^ v = 0, (2) 

where is the flow field in the fluid. The fluid is perfect if 
the energy-momentum tensor is diagonal in the local rest 
frame, i.e. viscosity and heat conduction are negligible. 
This can be assured if is chosen as 

= (e+p)u»u v -pg» v , (3) 

where e is energy density, p is pressure and g^ v is the met- 
ric tensor, diag(l, -1,-1,-1). The conservation equations are 



closed by the equation of state, which gives the relation- 
ship between e and p. Typically e = up is chosen, where 
the proportionality "constant" k may depend on temper- 
ature T, which is connected to the density n and pressure 
p via p = nT . In some solutions (as also on the analyzed 
one) a bag constant B can be introduced (this is not fa- 
vored by the data, because no first order phase transition 
is seen in high energy heavy ion collisions, these are in the 
cross-over regime of the QCD phase-diagram [3]). The ex- 
act, analytic result for hydrodynamic solutions is, that the 
hadronic observables do not depend on the initial state or 
the dynamical equations separately, just through the final 
state [1]. Thus if we fix the final state from the data, the 
equation of state can be anything that is compatible with 
the particular solution. This is the framework of several 
hydro solutions as detailed in the next paragraph. 

Even though many solve the above equations numeri- 
cally, there are only a few exact solutions for these equa- 
tions. One (and historically the first) is the implicit so- 
lution discovered more than 50 years ago by Landau and 
Khalatnikov [SHI]- This is a 1+1 dimensional solution, 
and has realistic properties: it describes a 1+1 dimen- 
sional expansion, does not lack acceleration and predicts 
an approximately Gaussian rapidity distribution. Another 
renowned solution of relativistic hydrodynamics was found 
by Hwa and Bjorken solution [8H10|: it is simple, 1+1 di- 
mensional, explicit and exact, but accelerationless. It is 
boost-invariant in its original form. Boost invariance is 
however incompatible with data from RHIC, so the so- 
lution fails to describe the data (it still can be used to 
estimate the energy density in high energy heavy ion col- 
lisions; note that this solution rather underestimating the 
initial energy density, see ref [1 lj ) . Important are solu- 
tions |111[T2"] which are explicit and describe a relativistic 
acceleration, i.e. combine the properties of the Landau- 
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Khalatnikow and the Hwa-Bjorken solutions. With these 
one can have an advanced estimate on the energy density, 
but investigation of transverse dynamics is not possible by 
these solutions. 

There were no 1+3 dimensional relativistic exact para- 
metric solutions investigated yet: the only exact solution 
is the one in ref. [2]. Observables from this solution were 
not computed and compared to data yet. Present paper 
hence calculates observables from a realistic 1+3 dimen- 
sional solution and compares them to data for the first 
time. There were several numerical solutions of relativis- 
tic hydrodynamics compared to data, but our method is 
different from numerical calculations: here one can deter- 
mine the best values of the parameters of the solution by 
fitting the analytic results to data. 



3 The analyzed solution 

The analyzed solution [5] assumes self-similarity and ellip- 
soidal symmetry. The ellipsoidal symmetry means that at 
a given proper time the thermodynamical quantities are 
constant on the surface of expanding ellipsoids. The ellip- 
soids are given by constant values of the scale variable s, 
defined as 



x(ty v{ty z(t) 



2 ' 



(4) 



where X(t), Y(i), and Z(t) are time dependent scale pa- 
rameters (axes of the s = 1 ellipsoid), only depending on 
the time t. Spatial coordinates are r x , r y , and r z . The 
velocity-field is described by a Hubble-type expansion: 



u^(x) =7 1 



X(t) Y(t) Z(t) r 
X(t) x, Y(t) y ' Z(t)' 



(5) 



where x means the four-vector (t,r x ,r y ,r z ), and X(t) = 
dX(t)/dt, similarly for Y and Z. The X(t) = X , Y(t) = 
Yo, Z(t) = Zq (i.e. all are constant) criteria must be ful- 
filled, ie. the solution is accelerationless. This is one of the 
drawbacks of this solution. 

The temperature T(x) and number density n(x) are: 



n(x) = n f— ) u(s), 
(*)=T p) 3/K * 

V T / V{s) 
, \ /T \ 3 ( K + 1 )/ K 

P(x) = Po {—) 



T 



(0) 
(7) 
(8) 



where r is the proper time, s is the above scaling vari- 
able, u(s) is an arbitrary function, while n$ = n| s= o jr=To , 
T = T| s=0 , r=To and p = p| s=0ir=To with p = n T 
(hence p does not depend on the spatial coordinates only 
r). Furthermore, ro is an arbitrary proper time, but prac- 
tically we choose it to be the time of the freeze-out, thus 
Tq is the central freeze-out temperature. Note that in our 



solution, the parameter n is arbitrary, i.e. any value of k 
yields a solution. The function v(s) is chosen as: 



-bs/2 



(9) 



where b = ^-\ r is the temperature gradient. If the fireball 
is the hottest in the center, then b < 0. 



4 Freeze-out and source function 

The picture widely used in hydro models is that the pre 
freeze-out (FO) medium is described by hydrodynamics, 
and the post FO medium is that of observed hadrons. 
Note that cases have been analyzed where the pre and 
post FO physical parameters are different, see refs. [T3T - H5] 
for details. In our framework we assume however that the 
freeze-out can happen at any proper time, e.g in case of 
a self-quenching effect or if the phase space evolution is 
that of a collisionless gas. Thus there is no jump in the 
equation of state post and pre FO, i.e. k goes to Kf roc 
smoothly, to the EoS of free hadrons. This is a widely 
used assumption, the framework of our solution is similar 

to that of refs. [THHM17]. 

In this case the hadronic observables can be extracted 
from the solution via the phase-space distribution at the 
FO. This will correspond to the hadronic final state or 
source distribution S(x,p). 

As mentioned in the previous section, one does not 
need to fix a special equation of state, because the same 
final state can be achieved with different equations of state 
or initial conditions [3]. Thus in this paper k is arbitrary 
- the hadronic observables do not restrict the value of k. 
The bag constant B does not have to be specified either, 
but the instantaneous FO we assume corresponds to zero 
bag constant. 

In our solution the source distribution takes the fol- 
lowing form: 



S(x,p)d 4 x = AfB(x,p)H(T)dTd 3 E^Pr, 



(10) 



where Af — g/(2-7r) 3 (with g being the degeneracy factor), 
H(t) is the proper-time probability distribution of the 
FO, B(x,p) is the Boltzmann-distribution and d 3 S^(x)p ti 
is the Cooper-Frye factor [18] describing the flux of the 
particles, and d^S^x) is the vector-measure of the FO 
hyper-surface. We assume that the FO happens at a con- 
stant proper time ro, i.e. H{r)dT = S Tq (t) and d 3 S^(x) = 
u fJ, d 3 x/u°, i.e. the FO hyper surface is assumed to be nor- 
mal to u^ L . 

We assume for the Jiittner-distribution: 



B(x,p) — exp 



/i(x) -p ll u^(x) 



T(x) 



n(x) exp 



PnU^ix) 



T(x) 



(11) 



where /i(x)/T(x) = hin(x) +^o/?o is the fugacity factor. 
Finally the source distribution is: 



S(x,p)d i x = Nn(x) exp 



p^u fi (x) 



T(x) 



P^ 



-S To (T)dTd 3 x. 

(12) 
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Note that the source distribution is normalized such 
as / S(x,p)d xd 3 p/E = N, i.e. one gets the total number 
of particles N (using c=l, h=l units) 



5 Observables from the solution 



Now let us calculate the observables that are usually mea- 
sured in high energy heavy ion collisions. The invariant 
momentum distribution Ni(p) can be calculated as: 



For the integration, we did a second order Gaussian 
approximation. The result for the invariant momentum 
distribution of a particle with mass m is: 



JVi(p) = / S(x,p)d 4 x = 2nNEV 
E 2 + m 2 



(14) 



x exp 



P 2 v 



2ET 2ET X 2ET y 2ET Z 

with the following auxiliary quantities: 

/o^r ^2^\ 3 / 2 
N = NnJ 2T ^ 



E=\E- 



V E 



(15) 



E 



E 



E 



V 



\-^\ d-S) 



Tr 



T,, 



Furthermore, T x ,T y , T z are the effective temperatures, i.e. 
inverse logarithmic slopes of the distribution: 



T x — To 



T z =Tn 



b(T - E) ' 
ET Y 2 

b(T - E) : 
ETpZ 2 

b(T — E) 1 



(18) 
(19) 
(20) 



where Xo, Yq and Zq are the (constant) expansion rates 
of the fireball, To its central temperature at FO and b 
the temperature gradient. There is an important criterion 



5.1 The transverse momentum distribution and the 
elliptic flow 

The invariant momentum distribution depends on all three 
momentum coordinates. Here we take the longitudinal mo- 
mentum to be zero, as the data we compare it to are mea- 
sured at mid-rapidity. We use transverse polar coordinates 
4> and pt instead of p x = p t cos(0) and p y = p t sm(<f>). This 
way Ni (p) can be rewritten as 



N 1 (p) = N 1 (p t 



v n cos(n0) 



(21) 



where v n are the flow coefficients, in particular vi is the 
elliptic flow. This way we may calculate the transverse 
momentum distribution Ni(p t ) from N±(p) of eq. (|14[) : 



and the elliptic flow: 



«a(Pt 



J£_ d0Ni(p) cos(20) 



(22) 



(23) 



Result for the transverse momentum spectrum is: (we 
substitute E = m t = y m 2 + p\ as being at mid-rapidity) : 



m(pt) V mt 



x exp 



P 2 (T eff -r ) 

m t T eS 
n z pf 



(24) 



2m t T 



2m t T r 



tJ-cfi 



where we introduced 1/T cff = 0.5(1/T X + 1/T y ), the effec- 
tive temperature. The result for the elliptic flow is: 



V2(Pt) 



h(w) 



(25) 



where Iq, and I\ are the modified Bessel functions while 



PI 



imt \T, 



1 



1 



(26) 



The formula for vi gives back previously found formulas of 
non-relativistic solutions [TB] and relativistic solutions |191 
[21?] . Also the formula for Ni(p t ) is similar to results of the 
previously mentioned papers. 

An important consequence of the above results is that 
neither Ni nor V2 depends on the EoS itself, only through 
the final state parameters. If we determine for example To, 



for the validity of this calculation: T x y x > Tq has to be the freeze-out central (at the center means here r x 



true. This is the case if b < (i.e. the fireball is cooler 
on the outside than in the inside) and E > Tq. In case 
of roughly 200 MeV central freeze-out temperature, this 
yields p > lAQMeV/c. We will see later, that this is true 
for all the data points we use. 



r, = 0) temperature, n or the initial temperature T; n iti a i 
still cannot be calculated. We only know that T; n i t i a i = 
To(To/T initia i) 3 / K , see eq. ©. Thus k or Initial has to be 
determined from another measurement, e.g. the spectrum 
of thermal photons. 



4 
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5.2 The Bose-Einstein correlations 



The two-particle Bose-Einstein correlation function [3T| is 
also calculable from the source distribution 



S(q,K) 



S{q = 0,K) 



(27) 



with q being the momentum difference and K the average 
momentum of the pair (and q <C K, hence p\ « p 2 ~ K), 
and S(q, K) is the Fourier transformed (x — > q) of S(x, K). 
The result of the solution is: 



C 2 (q, K) = 1 + cxp [-R 2 x q 2 x - R^q 2 - R 2 z q 2 z ] 



(28) 



where R x ,R y ,R z are the correlation radii, called HBT 
radii: 



R~, 



R l = 



Ri 



T t 2 (T x - 


To) 


M t T x 




T r$(T y - 


To) 


M t T y 




T t 2 (T z - 


To) 


M t T z 



(29) 
(30) 
(31) 



where M t is the transverse mass belonging to the average 
momentum K = 0.5(pi +P2), which is (at mid-rapidity) 
M t = 0.5 (m^i + TOt,2 )■ The m^i, and m t ^ quantities are 
the transverse masses, the T x . T y , and T z are the effec- 
tive temperatures belonging to the average momentum 
(i.e. here T x = T x \M t )- Note that the above formulas re- 
semble the usual scaling of the HBT radii: R 2 oc 1/M t . 
This means that all radii scale with the average trans- 
verse mass, independently of particle type, i.e. if plotted 
versus average mt, kaon radii and pion radii fall onto the 
same scaling curve. Note that HBT radii depend also only 
on final state parameters, not on the EoS itself, similarly 
to the elliptic flow or spectra (see the last paragraph of 
subsection 15. II on this matter). 

To compare the HBT radii with the data the Bertsch- 
Pratt [22] frame is to be used. It has three axes: the out is 
the direction of the average transverse momentum of the 
pair, the long direction is equal to the direction z, and the 
side direction is orthogonal to both of them. The result 
for i? out , i? s ide and i? long is: 



Ri 



-"side 



R„ 



R l 



n2 _ n2 
^Mong .-2 ' 



(32) 
(33) 



fit parameter 


Ni and HBT 


elliptic flow 




0-30% cent. 


0-92% cent. 


T [MeV] 


199±3 


204±7 


e 


0.80±0.02 


0.34±0.03 


u 2 t /b 


-0.84±0.08 


-0.34±0.01 


T [fm/c] 


7.7±0.1 




Zl/b 


-1.6±0.3 




NDF 


41 


34 


x 2 


171 


256 


X 1 with 3% 
theory error 


24 


66 



Table 1. The parameters obtained from fits to Au+Au 
PHENIX data |23H25| . The difference of the parameters can be 
explained by the centrality of the datasets (0-30% for spectra 
and HBT versus 0-92% for elliptic flow). The fact that u\jb is 
negative means that b < 0, i.e. the fireball is the hottest in the 
center and colder outside (i.e. a Gaussian temperature profile). 



6 Comparing the observables to RHIC data 

We compared the above results to PHENIX data of 200 
GeV Au+Au collisions. We fitted our above formulas to 
spectra and HBT positive pion data [231I2U (0-30% cen- 
trality) and elliptic flow data [25] for 7r ± , if ± , p and p 
particles (0-92% centrality). 

The used parameters for spectra and HBT are To (cen- 
tral freeze-out temperature), To (freeze-out proper-time), 
b (temperature gradient) and the expansion rates Xq, Yq 
and Zq. Instead of Xq and Yq however we use two more 
commonly used parameters, freeze-out expansion aniso- 
tropy e and average transverse expansion rate ut- 



X 2 ~Y 2 
X 2 + Y 2 ' 

-(— — 

2 V X 2 + Y 2 



(34) 
(35) 



. It turns out however, that the expansion rates appear 
only in combination with the temperature gradient b, hence 
we use u 2 /b and Z^/b as fit parameters. In eq. (f2~5l) . the 
formula for the elliptic flow, To and Zq do not appear, 
hence we do not use them when fitting vi data. 

For analyzing the confidence levels we calculate x 2 and 
number of degrees of freedom (NDF) of the fits. However, 
in the calculation there are numerous approximations. The 
main approximation is the Gaussian one in the integra- 
tion. Here we neglect terms of higher order. Later, with 
the given parameters we can estimate their contribution 
at a given point. The error coming from the integration 
can thus be estimated to be at least 3%. In table [T] we 
also give the x 2 using this additional error. 



Clearly in this solution the out and side radii are equal. 
This can be attributed to the instantaneous freeze-out; a 
non-zero freeze-out duration would make R 2 ut bigger by 
a term of Ar 2 p 2 /E 2 . Supported by the data, we use the 
At = approximation in our solution, which corresponds 
to instantaneous freeze-out. 



7 Discussion of the results 

The fit results are shown in fig. [T] The fit parameters are 
listed in table [T] The central freeze-out temperature To 
is around 200 MeV for both datasets, and the fireball is 
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colder away from the center. The expansion eccentricity e 
being positive tells us that the expansion is faster in-plane. 
Because of the Hubble-fiow this means that the source is 
in-plane elongated, similarly to the result of ref. [53] • The 
freeze-out happens at a proper-time of To = 7.7 fm/c. 

Our fit parameters describe the fireball at the freeze- 
out. However, the solution is time-dependent, most im- 
portantly the temperature depends on time as described 
by eq. ([7]) . We plotted the time-dependence of the central 
temperature in fig. [5] for several values of k, i.e. several 
EoS'. From this, assuming for example an average k of 
10 [27] one can also calculate the initial central tempera- 
ture of the fireball based on eq. 

/ \ 3 / K 

^initial = T I — I (36) 

\ ^initial / 




0.1 1 1 1 1 1 1 1 1 1 — 

200 400 600 800 1000 1200 1400 1600 1800 
p t [MeV/c] 



This yields 370 MeV at ti n itiai=l fm/c (note that t = r 
at the center of the fireball), which is in agreement with 
PHENIX measurements [25] . 

All the hadronic observables depend on the time of 
the freeze-out through the central temperature parameter 
Tq. If the freeze-out happens earlier or later, the central 
temperature is larger, and this enters into the HBT radii 
and the elliptic flow as well through the centeral freeze-out 
temperature. However, if we fix these final state (freeze- 
out) parameters and assume a k value, we can go back in 
time using the solution. We have to take eqs. ([25]) and (|32|) 
and substitute the time dependent central temperature: 



3/k 



(37) 



See such plots in fig. [2] for different values of k. In this 
figure we plot v 2 and R out = -Rsido for illustration purposes 
at pt = 400 MeV/c, the other parametes are taken from 
the fits. Note that in reality k changes over time, i.e. it 
goes smoothly to the EoS of a collisionless hadron gas. 



8 Summary 

Exact parametric solutions of perfect hydrodynamics were 
long searched for in order to describe the matter produced 
in heavy ion collisions at RHIC. We extracted observables 
for the first time from the relativistic, 1+3 dimensional, 
ellipsoidally symmetric, exact solution of ref. [2]. We cal- 
culated momentum distribution, elliptic flow and Bose- 
Einstein correlation radii from the solution. We compared 
the results to 200 GeV Au+Au PHENIX data [21H25] . The 
solution is compatible with the data. The fitted parame- 
ters of the solution describe the hadronic freeze-out. In 
the framework of our solution the fireball is in-plane elon- 
gated and has a Gaussian temperature profile. If using an 
experimentally determined average EoS of n w 10 [21], our 
results yield approximately 370 MeV at Ti n iti a i=l fm/c, in 
agreement with recent PHENIX measurements |28| . 




200 300 



400 500 600 
p t [MeV/c] 



700 800 




200 400 600 800 1000 1200 1400 1600 1800 
p t [MeV/c] 

Figure 1. Fits to 0-30% centrality PHENIX Au+Au spec- 
tra [23 (top) HBT radii [25] (middle) and 0-92% centrality 
PHENIX Au+Au elliptic flow [24] (bottom). See the obtained 
parameters in table [T] In the middle plot the lower curve is the 
fit to R ou t and R S ide- 
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Figure 2. Time dependence of the central temperature of the 
fireball, from eq. (J7J) (top) and time dependence of the HBT 
radii (middle) and the elliptic flow (bottom) where time depen- 
dence is included in the temperature To. The plots are shown 
for different k values. In reality k may change with time, we 
show here the curves only for fixed k values. Assuming an aver- 
age of k — 10 |27| onegets an inital temperature of 370 MeV at 
iinitiai=l fm/c, in agreement with PHENIX measurements |28| . 
The parameters used in this figure are that of table [T] Addi- 
tionally a transverse momentum pt had to be specified to plot 
V2, -Rout and -R a ide at, we chose pt = 400 MeV /c. Note further- 
more that as we go backwards in time, the central temperature 
reaches a point where the E > To criterion is not valid any 
more. Thus in V2, Rout and i? s idc we can go back only while 
T < E is valid, i.e. roughly 4 fm/c for k = 3 at the given pt- 
We did not plot the k = 1.5 case on those plots as the validity 
ends there at much bigger times. 
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